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Abstract 

Many problems give rise to polynomial systems. These systems often have several parameters and 
we are interested to study how the solutions vary when we change the values for the parameters. Using 
predictor-corrector methods we track the solution paths. A point along a solution path is critical 
when the Jacobian matrix is rank deficient. The simplest case of quadratic turning points is well 
understood, but these methods no longer work for general types of singularities. In order not to miss 
any singular solutions along a path we propose to monitor the determinant of the Jacobian matrix. 
We examine the operation range of deflation and relate the effectiveness of deflation to the winding 
number. Computational experiments on systems coming from different application fields are presented. 
2000 Mathematics Subject Classification. Primary 65H10. Secondary 14Q99, 68W30. 
Key words and phrases, deflation, Newton's method, path following, polynomial system, singular 
solution, sweeping homotopy. 

1 Introduction 

We consider systems /(x, A) = of polynomial equations / — (/i, /2, . . . , /at) in n variables x — 
{xi,X2, ■ ■ ■ ,Xn) and m parameters A = (Ai, A2, . . . , Am), with complex coefficients: fk € C[x, A], k = 
1,2, . . . , N. In this paper we restrict to isolated solutions so we assume N > n and most often N = n. 
For random choices of the parameters A, we expect all solutions to /(x. A) = to be isolated and well 
conditioned. Even though the total number of solutions may grow exponentially in the dimensions and 
degrees of the system, numerical homotopy continuation algorithms are efficient to enumerate all solutions. 

In many applications one wants to know for which values of the parameters A, the corresponding values 
of X satisfying /(x. A) = give rise to singular solutions. For values of A, the solution x is singular if the 
Jacobian matrix A(x), where Aij = for i = 1,2, . . . , N and j = 1,2, . . . ,n, has rank strictly less than n 
when evaluated at x. A value for A is critical if some corresponding values for x are singular solutions. To 
locate all critical values, we may solve an augmented system: 

/(x,A) = 

Fi^,X,ti)^{ A(x,A)m = (1) 

C^fl = 1 

where /x consists of n additional multiplier variables and c e C" is a tuple of n random complex numbers. 
The condition — 1 implies that we look for solutions (x. A, /i) to express as a nonzero linear 
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combination of the columns of the Jacobian matrix y4(x, A). The system in ([T]) is an example of the 
Jacobian criterion to find critical values. 

Solving F(x, A, /x) =0 is an effective global method to locate all critical values A, but because the 
input size of F is more than double the size of the original system /, this global method is often too 
expensive — its underlying complexity is that of the discriminant variety. In this paper we focus on the 
local approach: given a sufficiently generic value for A, we start at the corresponding solutions for x and 
trace the algebraic curves as we sweep through the parameter space. 

For introductions to homotopy continuation methods specific for polynomial systems, we recommend |18| . 
[24] , and [3^ . The books [I] , [9] , and [23 provide introductions to path following methods applied to general 
nonlinear systems and systems of differential equations. The authors of [7J study polynomial differential 
systems in the real plane and developed software to draw phase portraits. Computer algebra is used to 
compute all singularities, but it is noted in [7j that for high degrees this can take a long time. Recent 
related symboHc methods are described in [14] and [29] . 

In most applications, one is mainly interested in real solutions. However, a complex solution curve of 
a polynomial system may have isolated real solutions. Such a real solution on a complex curve will be 
isolated in the real space and will manifest itself as a singular solution on the curve. The methods in [52] 
rely on a global application of a deflation operator to locate real solutions on a curve in complex space 
while the sweeping homotopies defined in this paper offer a local approach for this problem. 

The contributions of this paper are twofold. To detect general types of singularities along a solution 
path we first describe an algorithm to monitor the determinant of the Jacobian matrix. Secondly, we 
investigate the effectiveness of deflation to accurately locate the detected isolated singular solution. These 
two contributions are outlined in sections 4 and 5. In sections 2 and 3 we first define sweeping homotopies 
and describe the problem statements. This paper ends with a report on our computational experiments. 

2 Sweeping Homotopies 

The system /(x. A) = defines already a homotopy. We call it a natural parameter homotopy because the 
parameters A appear naturally. To track the solution paths x(A) with pseudo arc length continuation we 
first compute a tangent vector v = (vx,va)"^ at the current point (xq, Aq) e C" x C": 

( ||(xo,Ao) f(xo,Ao) IMI = 1- (2) 

At a regular solution on an algebraic curve defined by a complete intersection where N = n and to = 1, 
the unit tangent vector is defined uniquely up to orientation. After selecting an orientation, a prediction 
for the next point is then (xi, Ai) = (xq, Aq) + /i(vx, vx) where h is the step size. Using interval methods 
as in [13] to control the step size, one can rigorously prevent jumping from one path to another. After the 
prediction step, Newton's method is applied to correct the predicted solution back to the solution curve. 
Algorithms for the adaptive use of multiprecision arithmetic during path following are proposed in [2] 
and [12]. 

In our sweep we use an artificial parameter homotopy, introducing a new artificial parameter t to make 
a convex combination between two given sets of parameter values Ao and Ai. Given /(x, A) = 0, Ao and 
Ai, a sweeping homotopy is defined as 

'^(-'^'*)-{ (l-t)(A-Ao)+i(/-^0 = O ^^[O'l]- 

For t = 0, we start a solutions of /(x, Aq) = and as t moves to 1, we sweep to the solutions of /(x, Ai) = 0. 
We use the same type of pseudo arc length continuation as above, with t as an added parameter, enforcing 
the orientation of the tangent vector so that t always strictly increases in value. 
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The distinction between natural and artificial parameter homotopies has profound consequences for the 
treatment of singularities. Consider for example f{x, X) — + — I — 0. Viewing f{x, A) = naturally, 
we recognize the equation of the (real) unit circle in the plane. Tracking the curve x(A) as defined by 
the natural homotopy f{x, A) = is then simply tracing the circle, either clockwise or counterclockwise. 
Consider picking Ag = and Ai = 2 as start and target values in a sweeping homotopy: 

'^(-'^'*)-{ (i-t)r+t(fr2;:2 (4) 

Forcing the orientation of the tangent vector to the path {x{t), X{t),t) so t is strictly increasing leads to 
a quadratic turning point for t = 0.5. At that point, the two real paths turn into the complex plane. 
While real homotopies (i.e.: homotopies with all coefficients real) for solving polynomial systems lead to 
singular points along the solution paths, as shown in [TH], generically, only a finite number of quadratic 
turning points occur. However, for our problem, perhaps we might assume a generic choice of the values 
for the parameters \q at the start of the sweep, but even that is insufficient to exclude general types of 
singularities. 

3 Detection and Location of Singularities 

Given a sweeping homotopy /i(x. A, i) = and a start solution (xq, Ao,0), our problem is then to detect 
and locate all singular solutions along the path as t advances from to 1. 

For the most common type of singularity, the quadratic turning points, the detection and location of 
singularities along the solution paths is done as follows: 

Detection: via the orientation of the tangent vectors The tangent vector v has three components 
V — (vx, va, Vt) and each time we force its orientation so vj > 0. If ti is before and ti after a turning 
point, forcing Vt > at ^2 will lead to a change in the angle between the corresponding tangent 
vectors v{ti) and v(t2) so that its inner product (v(ti), v(t2)) < 0. 

Location: shooting method for the step size Once we have two consecutive tangent vectors v(ti) 
and v(t2) for which (v(ti), v(t2)) < wc look to find h so that (v(ti), v(ti + h)) — 0. To find h we 
apply a shooting method. We found the description in |21| very clear and useful. 

The two tasks detection and location of a singularity along a solution path defined by a sweeping homotopy 
are more accurately described in the input / output statements in Table [1] and Table [2l 

Input/Output specification for Detection Problem 
Input : h(x, X,t) — sweeping homotopy, 
(xq, Aq, 0) a start solution. 
Output : solutions (x(ti), A(ti), ti) and (x(t2), A(t2), ^2), 
where the interval (ii,i2) contains a singularity. 

Table 1: The Detection Problem: detect singularities along a path. 

The output of the Detection Problem will be empty if there are no singularities for all t e [0, 1]. The 
main difficulty is to distinguish this case from the case where the solutions paths run straight through 
a multiple solution. If the paths are straight, the path tracker will not slow down and overshoot the 
singular solution. The output of the Detection Problem will be incomplete in the case where paths near 
a singular solution are very hard to follow. An example of that case is when a severe drop in the rank of 
the Jacobian matrix causes Newton's method to fail. The output in such a case then only consist of ti 
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and its corresponding solution, even as we have a good guess for t2, the solution corresponding to t2 will 
be missing. 

The difficulties in both problems are often complementary. On the one hand, singularities that are easy 
to detect (because Newton's method fails) , are usually hard to locate (again for the same reason) . On the 
other hand, multiple solutions for which Newton's method converges are hard to detect along a path. 

The output of the Detection Problem determines the input for the Location Problem, specified in 
Table H 

Input/Output specification for Location Problem 
Input : /i(x, A, i) = sweeping homotopy, 

(x(ii), A(ti), ti) and (x(t2), A(i2), ^2) are solutions 
where the interval (^1,^2) contains a singularity. 
Output : (x*. A*, t*) accurate singular solution of /i(x, A, t) = 0. 

Table 2: The Location Problem: locate accurately a singularity along a path. 

The quality of the input to the Location Problem will determine the difficulty of the Detection Problem. 
If the corresponding solutions at ti and t2 are accurate and not too far apart from each other, then solving 
the Location Problem will be much easier than when the corresponding solution for ti is inaccurate and 
the solution corresponding to t2 is missing. 

In the next two sections we outline our algorithms for the detection and location problem. 

4 Detecting Singularities along a Path 

The main difficulty in the detection problem is that the path tracker may not slow down when approaching 
a well conditioned multiple solution. Therefore, we first look for a criterion to decrease the step size and 
to back up towards previously computed values. 

The orientation of the tangent offers a clear criterion for quadratic turning points, but is no longer 
useful when sweeping paths that do not turn at a singularity. Monitoring the signs of the eigenvalues of the 
Jacobian matrix also captures many types of singular solutions, see e.g. jlOj . but we have encountered cases 
- see the applications section below - where the eigenvalues do not change signs when passing through 
the singular solution and where thus the determinant of the Jacobian matrix stays monotone of the same 
sign, only touching zero at the singular solutions. Our experiences have led us to opt for the determinant 
of the Jacobian matrix as the main criterion. This determinant is obtained as a relatively easy byproduct 
of the application of Newton's method. 

We can compute the determinant only we have as many equations N as unknowns n. In case N > n we 
can locally make the system square either by adding random combinations of the extra N ~n polynomials 
to first n polynomials, or by adding N — n slack variables to all polynomials, see [30l §13.5]. 

To detect singularities, we keep a window of three consecutive values for the artificial parameter t: 
ti < t2 < ts, along with the values of the determinants di, ^2, and of the Jacobian matrix at the 
corresponding solutions x(<i), x(t2), and x(i3). If there are any sign changes in the determinants, then 
the detection problem is solved. Otherwise, we compute an interpolating parabola p(i) so that p(tfc) = dk, 
for k = 1,2,3. If all determinants are positive, we compute the minimum of p. If all determinants are 
negative, we compute the maximum of p. The distinction on the sign of the determinant only matters for 
real solution paths. For solution paths in complex space, we monitor the modulus of the determinant, i.e.: 
dk = I det{A{x{tk), A(ife), tk)\, k — 1,2, 3. An exphcit criterion is given in ^ below. 
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Lemma 4.1 Let di, d2, and ds correspond to the three consecutive values for t: ti < t2 < t^. Then 



tjids - d2) + tl{di - ds) + ti(d2 - di) 

2{di{t2-t3)+d2{t3-ti)+d3{ti-t2)) ^' 

is a critical value for the interpolating parabola through the points (t2,d2) and (t^^ds). 

Proof. The Lagrange form of the parabola interpolating through the points (ti, di), (t2, ^2) and (ia, d^) is 

_ it-t2){t-t3) ^ (t-tl)ft-t3) , ^ {t-t,){t-t2) ^ 

^^'^ - {t,-t2){t,-t,)'^' + K^m^) ' ^ (ts-t,){t,-t2)'^'- 

Computing the derivative p'(t) = and solving for z in p'{z) — yields ([5]). □ 
Lemma 4.2 The cost of evaluating ([5]) along a path in C" is 0{n). 

Proof. Tracking a path, using Newton's method as a corrector, the Jacobian matrix is evaluated and 
decomposed when solving a linear system to obtain the next iteration of Newton's method. Given an LU 
decomposition of the Jacobian matrix, the extra cost of computing the determinant is just a product of n 
numbers along the diagonal of one of the factors L or U. Once the points {ti, di), (^2, (^2) Siiid [t^, ^3) are 
given, evaluating ([5]) requires a constant number of arithmetical operations. Thus the cost of evaluating ([5]) 
along a path in C" is 0{n). □ 



Lemmas 14.11 and 14.21 ensure that we have an explicit and efficient formula to evaluate along a path. 
Taking only the linear algebra operations into account, one iteration of Newton step costs at least 0{n^) 
which already dominates the 0(n) cost of evaluating the formula. To show the effectiveness of (O, we 
exploit the algebraic properties of our problem. In particular, we use Puiseux expansions ([B], [34j ) . 

Lemma 4.3 Assume the sweeping homotopy /i(x, A,t) = has an isolated finite singularity for t — t^. 
Then the determinant of the Jacobian matrix of A{:x.{t) , \{t)) fort sufficiently close to t^ equals (t — t^Y, 
for some fractional power p. 

Proof. A Puiseux expansion at an isolated singular solution of a polynomial system is a fractional power 
series, i.e.: with rational numbers for the powers in the series. In particular, for t ^ t^ we may write 
Xk{t) = ak{t — t^Y*'{l + 0{t)), k — 1,2, ... ,n, where ak & C \ {0} and Ofe e Q. Because we assumed the 
singularity to be finite (i.e.: not at infinity), we have that Ok > for all k. By definition of the sweeping 
homotopy, X(t) is linear in t and thus also in i — . 

Substituting Xk{t) = Ck{t — t^,)"-'' {1 + 0{t)) and the linear expression for Xk{t) into the Jacobian matrix 
A(x(<), A(t), i), we view A as A{t), as a matrix of polynomials with fractional powers in t — So also 
its determinant, det{A(t)) is a polynomial in t ~ t*. Ignoring higher order terms, we let p be the smallest 
power of i - in dct{A{t)). For t « t*, we then have det{A{t)) « (< - t*)^. □ 

If the power p of Lemma 14.31 is a natural number, then the interpolating parabola of Lemma |4. II will 
locally resemble very well the determinant itself and its critical value will be close to t*. For fractional 
powers of p, the determinant does not look like a polynomial. For example, consider t^ — and p = 1/2, 
then for t > 0, det(A(t)) = y/i and for t < 0, det{A{t)) = ^/^. 

Theorem 4.4 Assume the sweeping homotopy /i(x. A, <) = has exactly one isolated finite singular solu- 
tion at e [fijfa]. Then for any t2 G (tijt^) and z as in z G [ti,i3]. 
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Proof. Without loss of generality, consider t^, — 0. If p of Lemma is a natural number, then for i sa 0: 
det{A{t)) w tP and the interpolating parabola of Lemma[5]will have its critical value z in [tijia]. Even 
if p is not a natural number and a fraction, for all t2 G (ti, t^) the value of det{A{t2)) will be smaller than 
det{A{ti)) and det{A{t3)) so the interpolating parabola of Lemma[5]will have the right concavity and thus 
its critical value z will also be in [ti,t3]. □ 

Once the formula ([5]) yields a value z € [ti, t^] we then have a candidate singularity at t* and we 
need to locate it accurately. The case where the determinant of the Jacobian matrix has a local optimal 
value in [ti, ts] but no root is captured by the application of a minimization algorithm on the determinant, 
viewed as a function in t. Although we could further apply parabolic interpolation and use z of ([5|) as the 
next value for t, for fractional powers, we recommend the golden section search method [26]. The golden 
section search method requires an optimal number of function evaluation and is guaranteed to find the 
optimum if the function is unimodal over the interval. 

The key assumption of Theorem 14.41 is that there is only one singularity in the interval. For a reliable 
implementation of this algorithm, we must relate the step size h of the path tracker to the distance S between 
two singular solutions. If h is sufficiently smaller than d, then it is safe to assume that the determinant 
of the Jacobian matrix will behave as a unimodal function between three consecutive predictor-corrector 
steps. For polynomial systems, the total degree D is the product of the degrees of the polynomials in the 
system. Following Bezout's theorem it is a crude upper bound on the number of isolated solutions and 
therefore also on the total number of singularities. Assuming a uniform distribution of the singularities, a 
pessimistic lower bound on the step size h is 



5 Locating Singularities along a Path 

If a singular solution at t* is hard to detect, then for almost all t close to t* the Jacobian matrix is 
sufficiently well conditioned for Newton's method to converge well. In that sense, getting close enough to 
singularity to locate it with sufficient accuracy is then no problem. Thus then the main difficulty with the 
location problem occurs when Newton's method fails. 

The solution to the detection problem has made the location problem similar to an endgame |25j (see 
also [SI])- In this section we discuss the effectiveness of defiation — the idea to apply deflation to locate 
singular solutions of polynomial systems occurred first in [28] , see [15] for a symbolic deflation method — 
to locate general types of isolated singular solutions in the context of a sweeping homotopy. 

The deflation operator works locally starting at an approximation for a singular solution for which the 
Jacobian matrix has numerical rank is R. Numerical rank revealing algorithms can be found in '20J . Then 
R + I multiplier variables fi are used. To reduce to the corank one case, we multiply the Jacobian matrix 
with a random matrix B. This matrix B has R + 1 columns and as many rows as the columns of the 
Jacobian matrix. Then we apply Newton's method on the system 

r /(x,A) = 

E{x,X,fj,)=l A{x,\)Bti = (7) 
[ c^/x = 1 

The system E{x, X, fi) — looks very similar to the augmented system of the Jacobian criterion ([T]), with 
the addition of the numerical rank as extra local information. One application of the deflation operator 
may not be enough to completely recondition the isolated singularity and we have to apply deflation 
recursively. As proven in [5j and [16^, the number of deflations needed to restore the quadratic convergence 
of Newton's method is strictly less than the multiplicity of the singular solution. 

The term endgame operation range was coined in [25j . In general, this endgame operation range is the 
range for which the endgame techniques are effective. If we use for example extrapolation methods, then 
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we need on the one hand take samples along the path close enough to the singularity. On the other hand, 
if we get too close to the singularity, the iterations produced by Newton's method will be too inaccurate 
for the extrapolation. If we can adjust the working precision of our calculations, then we can guarantee 
that the endgame operation range is nonempty. 

The idea of deflation is to consider in addition of the original polynomials also the derivatives. If on 
the one hand, we are too far from the singularity, then adding the derivatives at the current approximation 
may lead to an inconsistent problem and lead to divergence. On the other hand, getting close enough to 
the singular solution may no longer be possible by the plain application of Newton's method. 

For deflation, one critical factor in its endgame operation range is the winding number. The winding 
number occurs as the denominator in the fractional power (or Puiseux) expansion of the solution path 
at the singular solution. The multiplicity of the solution bounds the winding number. The higher the 
winding number, the harder it could be for the derivatives to vanish in the proximity of the singularity. 
Proposition 15 . II formalizes the relationship between the winding number and the endgame operation range 
for deflation. 

Proposition 5.1 Let /i(x, A,t) — be a sweeping homotopy with an isolated finite singular solution for 
t — tsf with winding number uj. If for some component k: hk{x{t), X{t),t) is 0{t) for t « t*, then 
||i^/i(x(i), A(i),i) could in the worst case be 0{t^/'^). 

Proof. Without loss of generality we may assume that = 0, A(i*) = and x(t,) = 0. Expanding the 
solutions x(i) at = leads to fractional power series Xi{t) = Cii'"'/"(1 + 0{t)), iov i — 1, 2, . . . , n, where 
Ci S C \ {0} and Vi is a natural number. By definition of the sweeping homotopy, the relation between A 
and t is simply linear and it is straightforward to express X{t) as a linear function of t. 

Substituting A(i) and the power series for x(t) in hk we obtain /ifc(x(t), A(i), i) ~ 'ytP{l + 0{t)) for 
some nonzero complex constant 7 and some power p. Since we assumed a finite singular solution: p > 1. 
Similarly, for a derivative we obtain ^^(x(t), A(i), t) = 5t'^{l + 0{t)) for some nonzero complex 
constant 5 and some power q. Because not all derivatives will vanish at the singular solution, suppose k 
and i are such that q is the lowest positive exponent. Take then p = [lo + l)/uj and \ci q ~ p — 1. □ 

As a practical result of Proposition 15.11 we may make some pessimistic predictions on the endgame 
operation range of deflation. For example, ii lj — A and we need the residual of the derivatives to be about 
10^^, then the residual of the approximation must be about 10^^. 

Also the numerator of the exponent in the leading term of the Puiseux series plays an important 
role as it determines how sharp the curve bends as t approaches i,. Even with a high winding number, 
extrapolation methods will be effective for low numerators, but as the numerator is close to oj itself, then 
the solution curve will appear to be linear unless we get really close to i*. 

To estimate the winding number, Richardson extrapolation (see e.g.: [3]) cannot be applied directly 
because the exponents in the power series are unknown — we refer to [1] for general extrapolation methods 
for unknown exponents. In extrapolation methods to estimate the winding number for diverging 
solution paths were developed. An algorithm to predict the order of the deflation was recently presented 
in [17]. 

6 Computational Experiments 

We have implemented the detection criterion in the publicly available open source software PHCpack [33] 
and applied it to three polynomial systems, coming from different application fields and documented in 
the literature [8], [22], [^. 
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6.1 a system from molecular configurations 



The following system occurs in [8]: 



/(x,A) 



i(a;| + 4x2X3 +xl) - 


1- A(X2X3 


-1) 


= 


i(x§ +4X3X1 +Xl) ' 


K A(xix2 


-1) 


= 


|(Xi + 4xiX2 + X2) - 


1- A(xf X2 


-1) 


= 



(8) 



The system is listed as a nontrivial example in [351 pages 391-392]. However, the system is small enough 
for global methods. The Jacobian criterion ([1]) gives a system we solved with the blackbox solver of 
PHCpack 33J. The system has 54 generic solutions which can be divided into five groups with the same 
xi, X2, X3 and A values. The first four groups have the same absolute values of xi, X2 and X3 with the 
natural parameter A being either +1.5i or — 1.5i, i = T. There are exactly twelve solutions in each of 
the first four groups. The last group corresponds to the approximate value ± 0.866025403780023 as the 
natural parameter A. For these two values there are curves of degree six. In this example, all critical values 
for A were found via the Jacobian criterion. 

As A approaches zero, the system becomes singular. At the origin, there is one solution of multiplicity 8 
for the system when the deflation method in PHCpack is applied. To test our detection algorithm, we 
consider sweeping A through zero. The sweeping homotopy is 



/(x,A) 



i(xi + 4x2X3 + x§) + A(x|xi - 1) = 
i(xi + 4x3x1 + xj) + A(x§x? - 1) = 



i(Xi + 4xiX2 + X2) + A(xfx2 

(A-l)(l-i) + (A + l)t = 0. 



1) = 



(9) 



As the artificial parameter t goes from to 1, the natural parameter A is swept from +1 to — 1. 
According to the multihomogenous Bezout bound [301 , the permanent of the degree matrix of the system 
is 16 for nonzero values of A. This bound is sharp, so all solutions in a multihomogeneous homotopy 
converge to finite solutions. Among the 16 solutions, four are symmetrical complex conjugate solution 
pairs and four are symmetrical real solution pairs. By the symmetry, the solutions break up into orbits 
of type xi — X2, X2 = X3, xi = X3 and xi = X2 = X3. As A is swept from +1 to —1, starting with start 
solutions at t = 0, four real solution paths converge around the origin and the four complex solution paths 
diverge. If we would like to track the converging complex solutions paths, we could set the homotopy to 
(A + 1)(1 — + (A — l)t = such that the four complex solution paths converge around the origin and the 
four real solution paths diverge. The special value zero for the natural parameter A is found by the sweep 
as the tangent flips. A solution of multiplicity 8 is found at the origin. 

Our new detection algorithm is needed to detect the singularities at A = ±0.866025403780023 for which 
there are curves of degree 6. Because even close to this critical value, the solutions are still relatively well 
conditioned, monitoring just the orientation of the tangent is insufficient. 



6.2 modeling neural networks 

Families of polynomial systems often not only depend on parameters, but also the dimension may scale. 
Our next example originates from [27] . For n = 3, an example of a system in this family is 

{X1X2 + xix| — Axi + 1 = 
X2XI + X2X§ - AX2 + 1 = (10) 
X3X1 + X3X2 — Ax3 + 1 = 0. 
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The application of the Jacobian criterion in ([T]), results in a 7-by-7 system with 54 regular solutions. 
Critical values for A found among these solutions are 0, and the approximations 1.88988157484231, 
3.61703146124952, 2.38110157795230, and -0.414704714645147. As the dimension n grows, the plain 
application of the Jacobian criterion will quickly lead to an intractable problem, whereas the complexity 
of tracking one solution path scales much better. 

The singular solutions for the critical values corresponding to A = were the hardest to detect and 
stimulated the development of our detection algorithm. The homotopy which sweeps A through zero is 

X1X2 + xix\ — \xi + 1 = 

, , . X2XI+X2XI-\X2 + 1 = Q 1 . 

J^'^^^'-S X3XI+X3XI-\X3 + 1 = ^^^^ 

(A + 0.1)(l-t) + (A-0.1)< = 0. 

As the artificial parameter t goes from to 1, the natural parameter A is swept from —0.1 to 0.1. Passing 
through A = 0, the tangent vector does not flip back, the determinant does not change sign and comparing 
the signs of eigenvalues for A < and A > does not reveal anything. Without our detection algorithm, 
the path tracker will not back up and the solution of multiplicity four for A = would go undetected. For 
general values of n, the solution corresponding to A = has multiplicity n + 1. 



6.3 a symmetric Stewart-Gough platform 

A Stewart-Gough platform consists of two plates connected by six legs. One plate is fixed (the base plate) 
while the other plate (the top plate) moves as the leg lenghts change. These platforms are use used for 
example in flight simulators. At a singularity the trajectory of the top plate is no longer uniquely defined. 
In our experiments, we follow [35] where the equations for a symmetric Stewart-Gough platform are: 

{x, - Xm)"^ + {y, - y,o)^ + zf -If =0,i = 1, 2, 6 
f(^,\-J {X2- xi)^ + {y2-- yi)'^iz2- ~2Rl{l- cos{ai)) = , 



{xi - xoy + (yi - yoV + (^i - ^0)^ - i?f = 

(2^2 - Xo)^ + {y2 - yo)^ + {Z2 - Zo)^ - i?f = 



for 



where 




'wf'xi+W^'wf'X2 

I mo 1 mo mi /-i n\ 

w^a 2/1 + ^2 w'a y2 (13) 

mo I mo mi 



w\ — 



3 sin(ai ) + (-l)'"\/3(cos(ai)-l) 
2 sin(ai ) 



w3 = 



(-l)"'^/3 
2sin(Qi) 




2 = -sin(ai-(-ir73cos(ai) ) - - ^, ' - - (^4) 

2sin(ai) I — n-m, — 1 f^T- T — R V / 



The polynomial system has three fixed parameters: ai, a2, and i?i which determine the configuration 
of the platform. The angle ol\ is the relative angle between two the triangles connecting the joints in the 
moving top platform, while 012 is the relative angle between two triangles connecting the joints in the fixed 
base platform. The radius i?i is the radius of joints on the top plate. As in [35J, we fix the configuration 
parameters: i?i = 1, ai and a2 are respectively 28 and 22 degrees. Although we could consider the system 
as depending on six parameters, the leg lengths Z^, i — 1,2,. ..,6, for the purpose of simplicity, we only 
treat Zi as a natural parameter. 

The symmetrical platform gives rise to a system of nine polynomial equations in nine unknowns x = 
{xo,yo, zo,xi,yi, Zi,X2,y2, Z2)- In the application of the Jacobian criterion, we need to solve a 19-by-19 
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polynomial system. Fortunately the system is sparse and the mixed volume of the tuple of Newton 
polytopes equals 4,608. Tracking 4,608 paths yields 256 regular solutions. 

Applying our sweep to find critical values is certainly much less expensive for this system. Fixing li to 
1.5, 2.0, and 3.0, we found four special values for the natural parameter li for each 1^ with higher precision 
than what was reported in [3S]. In addition, we are able to see that Zq can be either posive or negative. 
When li is around 1.003, a multiple singular point occurs at the origin and li approximates the value of 
Zi, the system becomes a two-parameter problem and requires special care. 
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